function [model,momentStru,KFStru,specStru,counterOut,names]=...
    modelAnalysis(modelHandle,param,model,sampleStru,reportStru,counterStru)
% =========================================================================
% modelHandle:  handle to a function 
%
% param:        parameter vector 
%
% model:        structure 
%
%       model.solveopt 
%       model.addsol 
%       model.data 
%       model.trainvec 
%
% sampleStru    structure 
% 
% reportStru 
% 
% counterStru   innovations  
%               aZero 
% =========================================================================

%% compareSrtru: Assign default values for the moments 
[~,reportStru]=ch_field(reportStru,'NSimVec',[100 size(model.data,1)]);
[~,reportStru]=ch_field(reportStru,'NRep',100);
[~,reportStru]=ch_field(reportStru,'NCorrel',4);
[~,reportStru]=ch_field(reportStru,'NIRF',8);
[~,reportStru]=ch_field(reportStru,'NVarDec',0);
[~,reportStru]=ch_field(reportStru,'flagUnitIRF',0);

%% Solve model 
disp(' '); 
disp('_______________________________________________'); 
disp('Solving Model'); 

[model.GG,model.RR,model.CONS,model.eu,model.cholV,model.ZZ,~,model.steady,~...
    ,names.steady,names.states,names.shocks]...
    =feval(modelHandle,param,model.solveOptions,model.addsol);

%% reportPositions: obtain positions of all objects to report in reportStru
[reportStru,reportFlag]=reportPositions(reportStru,names,sampleStru);

%% momentBase: Compute all moments 
momentStru=momentsDraws(model.GG(:,:,end),...
    model.CONS(:,end),...
    model.RR(:,:,end),....
    model.cholV(:,:,end),...
    model.ZZ(:,:,end),reportStru);

%% Counterfactual
if nargin > 5 && isempty(counterStru)==false
    disp('Performing Counterfactual');
    if isfield(reportStru,'counter')==1 && ...
            isempty(reportStru.counter.lessShocks)==false
        disp('Deleting shocks');        
        tempPos=cellposition( reportStru.counter.lessShocks, names.shocks.long );
        counterStru.innovations(:,tempPos)=0;        
    end
    
    [counterOut.yobs,statesCounter]=feval(model.addsol.counter,param,modelHandle,model.data,...
        model.solveOptions,model.addsol,...
        counterStru.innovations,counterStru.aZero);
    counterOut.states=statesCounter(:,reportStru.stateDecompPos);
    %counterOut.states=KFStru.smoothSt(:,reportStru.stateDecompPos);
else
    counterOut=[];
end

%% Smoother 
if isempty(counterOut)==true && isempty(reportStru.stateDecomp)==false;
    disp('Begin smoother');
    [KFStru,logL]=feval(model.addsol.funcKfilter,param,modelHandle,model.data,model.trainVec,...
        model.solveOptions,model.addsol);
    KFStru.analysis.states=KFStru.smoothSt(:,reportStru.stateDecompPos);
    KFStru.analysis.stateDecomp=KFStru.countSt(:,reportStru.stateDecompPos,:);
else
    KFStru=[];
end

if isempty(reportStru.stateSpectrum)==false;
    disp('Begin Spectrum'); 
    [specStru.cell,specStru.table,specStru.xlsCell]=estimaSpectrum(modelHandle,param,model.solveOptions,model.addsol,...
        reportStru.stateSpectrum,sampleStru,1000);
else 
    specStru=[]; 
end 